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The resonant state expansion (RSE), a novel perturbation theory of Brillouin-Wigner type de- 
veloped in electrodynamics [Muljarov, Langbein, and Zimmermann, Europhys. Lett., 92, 50010 
(2010)], is applied to planar, effectively one-dimensional optical systems, such as layered dielectric 
slabs and Bragg reflector microcavities. It is demonstrated that the RSE converges with a power 
law in the basis size. Algorithms for error estimation and their reduction by extrapolation are pre- 
sented and evaluated. Complex eigenfrequencies, electro-magnetic fields, and the Green's function 
of a selection of optical systems are calculated, as well as the observable transmission spectra. In 
particular we find that for a Bragg-mirror microcavity, which has sharp resonances in the spec- 
trum, the transmission calculated using the resonant state expansion reproduces the result of the 
transfer/scattering matrix method. 



PACS numbers: 03.50.De, 42.25.-p, 03.65.Nk 



I. INTRODUCTION 

Recently, a novel perturbation method for the treat- 
ment of open electromagnetic systems, the resonant 
states expansion (RSE), has been formulated. 1 Unlike 
previous perturbative approaches^— which due to their 
complexity and poor convergency are limited to small 
perturbations, this method is shown to be suitable for 
perturbations of arbitrary strength and shape. It is based 
on the concept of resonant states (RS) of an open system, 
also known in quantum mechanics as Gamow^ or Siegert 
states4i£ These states exponentially decay in time and 
grow in space at large distances, 10 Owing to their com- 
pleteness inside a finite area of space, RS can be used for 
expanding solutions of the Maxwell equations, reducing 
the wave problem for the modes of the system to diago- 
nalization of a finite complex matrix. Hence the strength 
and the shape of the perturbed dielectric profile that can 
be treated is limited only by the size of the chosen finite 
basis of unperturbed RS. 

The idea of resonances is a cornerstone of physics, al- 
lowing to rationalize the dynamic behavior of physical 
systems. In open systems, excitations decay with time, 
endowing resonances with a spectral width. Depending 
on their width and separation resonances appear in mea- 
sured spectra as isolated lines or merge into a contin- 
uum. The concept of RS provides a unified picture of 
an open system which includes all types of resonances 
and is an alternative to the commonplace division of the 
spectrum into non-decaying bound and continuum states 
with real energies. RS are discrete eigenstates which 
have complex frequencies (and equivalently energies and 
wave numbers) and satisfy outgoing wave boundary con- 
ditions. This corresponds to a physical situation that an 
open system, excited at an earlier time, looses its energy 
to the outside space. The imaginary part of the frequency 
reflects the temporal decay of the energy in the system. 
Owing to this leakage, the RS wave functions have tails 
(outgoing waves) which grow exponentially outside of the 
system and cannot be normalized by the usual integra- 



tion of their square modulus. Instead, the normalization 
and orthogonality of RS is given by an integral over the 
finite volume of the system and the energy flux to the 
outside in the form of a surface term; 



The presence of a continuum in the spectrum of a sys- 
tem is a significant problem for any perturbation the- 
ory. In open electromagnetic systems such a continuum 
is often the dominating if not the only part of the spec- 
trum. However, going away from the real axis to the 
complex frequency plane, the continuum can in many 
cases be effectively replaced by a countable number of 
discrete RS which form a complete basis. Therefore, RS 
of a perturbed system can be expanded into the unper- 
turbed RS. The expansion coefficients can be found by 
diagonalizing a complex symmetric matrix which con- 
sists of a diagonal matrix representing the bare spectrum, 
and the perturbation^ The perturbed resonant states 
can then be used to calculate the Green's function of 
the system via its spectral representation ; 12 i 13 using the 
Mittag-Lefffer theorem. The Green's function provides 
the complete system response and allows to calculate ob- 
servables such as emission, scattering, or transmission. 
This recently formulated general method of reducing the 
Maxwell equations for an open optical system to a lin- 
ear matrix eigenvalue problem is called resonant state 
expansion (RSE). 

The RSE has been suggested^ as an appropriate tool 
for calculation of sharp resonances in optical spectra, 
such as perturbed whispering gallery modes of a dielec- 
tric microsphere. Popular computational techniques in 
electrodynamics, such as the finite difference in time 
domain (FDTD)ii^ or the finite element methodic 
adapted to such problems, require a large computational 
domain in time and/or space, and can produce spurious 
solutions^ In particular, sharp resonances are charac- 
terized by optical modes which decay slowly in time and 
hence FDTD needs a large time domain. Furthermore, 
being applied to open systems, the finite element method 
either introduces a significant error when the boundary is 
too close, or needs to consider an excessively large domain 
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in real space in order to describe the far-field asymptotics 
correctly. The RSE does not suffer from these problems 
because it produces the eigenstates of the system, and in 
particular their wave numbers, directly by diagonaliza- 
tion of a matrix determined by the near-held properties 
only. 

In this paper, we apply the RSE method outlined 
in Section [TT] to various planar optical systems. For 
such effectively one-dimensional (ID) systems efficient 
alternative methods exist to calculate the transmission 
and reflection, enabling verification of the RSE results. 
We investigate the accuracy with which eigenfrequencies, 
eigenfunctions, the Green's function, and transmission 
can be reproduced. We give a method to evaluate the 
convergence and to extrapolate the results in Section ITlTl 
We apply the RSE to a perturbed dielectric slab in Sec- 
tion |IV] with two different kind of perturbations: a wide- 
layer perturbation, and a (^-perturbation. We find that 
the RSE converges to the exact solution with a power 
law in the basis size. As an example of a structure with 
a sharp resonance, we treat a Bragg mirror microcavity 
in Section MM 



II. THE PERTURBATION METHOD 

The system of Maxwell's equations for a planar dielec- 
tric structure with permeability /i = 1 surrounded by 
vacuum is reduced to the following wave equation: 

d 2 z E u (z,t) = [e(z) + Ae(z)] d t 2 E, y (M) , (1) 

where e(z) denotes the unperturbed dielectric profile, and 
Ae(z) the perturbation of the dielectric function. Here 
the transverse eigenmodes with index v are taken with 
zero in-plane wave number. The electric field E^(z, t) can 
be written in an harmonic form 



£j,(z,t) = £ v (z) exp(— icxyt) 



(2) 



with complex frequency cx v (c is the speed of light in vac- 
uum) and amplitude £ v (z) satisfying the following time- 
independent wave equation: 



{d 2 z + [s(z) + As(z)]4]£„(z) 



= 0. 



(3) 



The electric field £ u {z) and its first derivative arc con- 
tinuous everywhere. RS are eigenmodes which satisfy 
outgoing boundary conditions, given by the form 



£ v {z) = A v exp(ix„|z|) 



(4) 



in the surrounding vacuum with the amplitudes A~ (on 
the left-hand side) and (on the right-hand side of the 
structure) which are generally different. In the case of a 
mirror-symmetric system, A~ — for symmetric and 
A~ = — for antisymmetric modes. 

In the following, for the unperturbed RS, i.e. for 
Ae(z) = 0, £ v (z) is denoted as E n (z), and >c v as k n . 



The unperturbed RS are orthogonal and normalized ac- 
cording to 



e{z)E n (z)E m (z) dz 
E n (-a)E m (-a) + E n (a)E m (a) 



= S m 



(5) 



where z = ±a are the positions of the boundaries of the 
unperturbed system. The perturbed states are written as 
linear combinations of the normalized unperturbed RS, 



E n (z) 



(6) 



resulting in the linear eigenvalue problem for n v and c Ti 



v„ 



with the perturbation matrix 1 

V nm = / Ae{z)E n (z)E m (z)dz . 

J —a 



(7) 



(8) 



In a dielectric system with real refractive index, the RS 
wave numbers k n have the following general property: 
Imfc n < and Refc_„ = Refc+„. Additionally, in fD 
systems (planar systems at normal incidence) there is 
always a RS with Re fco = and Im ko ^ 0. We number 
the RS with increasing real part of their wave number, 
numbering the state with zero real part as state number 
zero. The number of RS in the unperturbed or perturbed 
systems is countable infinite. Therefore we always deal 
will a truncation of the basis of the RS, which is the only 
approximation of the theory. We refer to n max as the 
truncation number for the basis so that — n max < n < 
n mH . Hence the basis size N is given by 



N = 2n T 



1 . 



(9) 



Ideally, by choosing the basis size TV sufficiently large, 
the results of the perturbation theory can be produced 
with any given accuracy. 

The unperturbed system can be any convenient sys- 
tem. In the discussed ID case, a dielectric slab in vacuum 
having thickness 2a and real dielectric constant 



e(z) 



e s for \z\ < a 
1 otherwise 



(10) 



is the simplest system having an analytic solution. We 
use it as unperturbed system in the following. The ex- 
pressions for the unperturbed RS are given in the Ap- 
pendix. The dielectric constant is taken to be ^fk~ s = 1.5 
unless otherwise stated. 
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III. CONVERGENCE AND EXTRAPOLATION 

We introduce a method to estimate the convergence 
and to extrapolate the RS wave numbers calculated 
via the perturbation theory xi N ^ to their exact values 
^_(cxact) ^ we approximate the absolute error in 

each wave number as a power law in the basis size N: 



(exact) _ (AT) _ K ' N K „ K"N< . 



(11) 



We assume that the exponent in the power law (a' v or a") 
is a real number, so that the RS wave numbers converge 
in a straight line in the complex plane. To determine the 
coefficients and exponents of the two representations in 
Eq. (|f f I) we use four different values of N: Ni < N2 < 
N3 < N4, where 



N 2 a V 2 N 4 , iV 3 « r]N 4 , 



(12) 



producing four sets of wave numbers, with the factor 
< r\ < 1. We match states between the four sets se- 
quentially, i.e. first {xi N4 ^} to {xi N ^}, then {xi N ^} to 

{xi N2 ^}, and finally {xi^} to {xi }. In doing this, 
we use the following matching algorithm (MA) between 

two sets of wave numbers, {xi A ^} and {x^}: 

(a) Determine the distance between the complex wave 
numbers of all pairs with one element from {xi A ^} 
and one element from {x^}. 

(b) Select the pair with the shortest distance, store it, 
and remove it from the sets. 

(c) Repeat (b) until {xl^} or {x^} is empty. 



This procedure 

(Ni) (N 2 ) (Ns) 



results 



in N% vectors 
(xl lyi >,xl" 2, ,x^ 3, ,x^ i >) of RS wave numbers. 
The specific factors chosen between Ni, AT 2 , N 3 , and N4 
allow for the following analytical expressions for two sets 
of coefficients and exponents in Eq. (fTTj) . for each state 
v: 



In order to test the quality of our power law fit, we 
estimate for each state v the relative extrapolation error 
defined as 



F„=$(K' V N?»,KZN2-), 
where 2$(X, Y) = T(X, Y) + Y{Y, X) and 



T(X, Y) 



X 

1 

Y 



(18) 



(19) 



Indeed, F v has the meaning of a relative error in the 
power law approximation of the distance ^ oxact ' — xl N ^ 
deduced from the two sets of power law parameters. If 
this error is sufficiently small, F v < F max , and the power 
law converges (ot v < a max ), we can improve the result 
calculated for the largest basis size N4 by extrapolating 



it towards the exact value, xii 



(JV 4 ) 



—> x, 



(00) 



where the 



extrapolated wave vector x v is defined according to 
Eq. (JTTJ> as 



K,,N? 



(20) 



Otherwise, the power law is not describing the conver- 
gence well. We then use the absolute variation scaled to 
the system size to evaluate if the state has sufficiently 
converged 



Mi/ = max 

i=l,2,3 



ANa) 



A N >. 



(21) 



We use state v for the calculation of the Greens func- 
tion if its relative or absolute error is sufficiently small, 
i.e. if one of the two selection criteria (SC) is met: 

1. extrapolation error F v \K v N2"\a < M max provided 
that F v < F max and ct u < a max ; 

2. absolute error M v < M max . 

For the results shown in the present paper we used 



1, a„ 



-0.5, and r) = 2" 1 / 4 . 



K'„ 



K" = 




(13) 
(14) 
(15) 
(16) 



For extrapolation of eigenvalues and estimation of errors 
we also introduce mean values a„ and K v defined as 



a. 



K\,x a A" + KK 



(17) 



IV. RESULTS 

A. Wide-layer perturbation 

The perturbation being considered in this section is 
jiven by 



Ae(z) = 



Ae for a/2 < z < a , 
otherwise 



(22) 



with Ae — 10. The profiles of the unperturbed and per- 
turbed dielectric constants are shown in Fig.[T] The an- 
alytic solutions of the time-independent Maxwell's equa- 
tions using the RS boundary conditions are given in Ap- 
pendix El both for the unperturbed and the perturbed 
systems, along with the matrix elements V nm of the per- 
turbation. 
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FIG. 1: Dielectric constants of the unperturbed slab e{z) and 
a slab with a wide perturbation e(z) + Ae(z). The distance z 
is in units of the half width a of the slab. 



Using the procedure introduced in Section IIIII we cal- 
culate four sets of perturbed wave numbers and extrapo- 
late >r„ according to Eq. (|20|) . We also calculate the exact 



wave numbers xl° xact ' and match up exact and perturbed 
states using the MA. The resulting exact and extrapo- 
lated eigenvalues *c£°^ are shown in the inset of Fig. [21 
together with the unperturbed wave vectors. We measure 



the errors in xi relative to 



(exact) 



and compare it with T{xi Ni \ xr£, cxact ' 



byT^V^ ') 
to evaluate the ex- 
trapolation method. The results are shown in Fig.[5J We 
see that the relative error of the RS wave number is gen- 
erally reduced by extrapolation by more than one order 
of magnitude. 
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FIG. 

r(x& 

the RSE for the perturbation shown in Fig.[T] with and 
without extrapolation, respectively, for iV4 = 801. Inset: 
unperturbed and perturbed RS wave numbers; the latter 
are calculated analytically (empty squares) and via the RSE 
with extrapolation (crosses). 
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FIG. 3: Power law parameters and error estimates for the 
wide perturbation. (a),(b): Relative extrapolation error F v . 
(c),(d): exponent a v in the power law fit. (e),(f): absolute 
errors K v N a " and (g),(h): M v as function of the the state 
number u, calculated for different basis size N. The right 
panels display the data versus the state number v normalized 
to its maximum value n max = (N — l)/2. Straight magenta 
lines are a = —3 (c),(d) and power law fits (h),(f). 



The coefficients and exponents of the power law fit 
give us information about the convergence properties of 
the perturbed RS. For the wide perturbed layer they are 
shown in Fig. [3] We see in Fig.[3](a) that states close 
to the origin in complex wave number space (and hav- 
ing small state number values) are not described well 
by the power law (F u is larger than F maK ), even though 
Fig. [5] suggests that these states are well converged. This 
is reflected in the small absolute error M v shown in 
Fig.EIg),(h), passing the SC. We also see that for higher 
wave-number states passing the relative SC the exponent 
in the power law is close to a = — 3 [horizontal lines in 
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FIG. 4: (a) Wave number of the perturbed state v = 63 
calculated with different basis sizes N and extrapolated to 
the exact value, (b) Absolute "exact" error \x£^ c ^ — 
for different N (squares) and a power law fit (dashed line). 



Fig.[3](c),(d)], in accordance with the findings in Ref.Qi 

Furthermore, the absolute errors K u N a " and M v show 
universal dependencies on the normalized state number 
i//n max , as shown in Fig.[3](f) and (h). This provides us 
with a scaling law of the absolute errors versus the state 
number: 



M v oc {y/Nf 



(23) 



This cubic scaling is shown in Fig.[3](f),(h) by straight 
magenta lines. The power law exponent a also shows 
a universal dependency on the normalized state num- 
ber, being a = —3 for u/n max < 0.2 as can be seen in 
Fig. 0(d). In this region the states pass the relative SC 
and are extrapolated. 

An example of how the power law is applied to extrap- 
olate the wave number of a particular state ^ = 63 is 
given in Fig.@](a). Clearly, the extrapolation leads to a 
considerable improvement of the accuracy compared to 
wave number calculated with the maximum matrix size 
A/4. This is due to the good power law convergence as 
shown in Figure [4](b), seen by the straight line connect- 
ing the "exact" errors |>fl° xact ' 1 — xi N ^\ for the four basis 
sizes. 

The exact errors are only available if the exact solution 
is known, but in this ideal case we do not need the RSE. 
In a realistic case for which no such solution is known, we 
need to estimate the error of the power law extrapolation, 
which we do using the extrapolation SC and Eq. (fT8|) . In 
order to check how good this estimation is, we compare 
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FIG. 5: "Exact" relative extrapolation error 7?( exact ) versus 
relative extrapolation error F v for both accepted and rejected 
states, for N4 = 80f . The blue dashed line shows the antici- 



pated behavior F„ 
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F v with the exact relative extrapolation error pf; exact ^ = 
AT"" , xrl° xact ' ) — xi Ni ^). Such a comparison is shown 
in Fig. [5] for all states with a v < —0.5. We can see that 
the exact error ^7 , j oxact ) j s typically overestimated by F v , 
and for all states with F v < F max we have ^ cxact ) < ^ 
i.e. the extrapolation is improving the error. F v can thus 
be used reliably to verify the convergency and power law 
extrapolation. 



B. Electric fields 

The electric fields (EF) £ v {z) of the perturbed RS cal- 
culated via the exact formula Eq. (|A5[) are shown in Fig. [5] 
for a few lowest states in comparison with E n (z), the EF 
of the unperturbed RS, given by Eq. (|A1[) . The perturbed 
RS are normalized as in Eq. ([5])- In particular, their or- 
thonormality condition reads 



e p (z)£ 1/ (z)£ fi (z) dz 
£„(-g)£ fl {-g) + £ u {a)£ l _ l {a) 



(24) 



where e p (z) = e(z) + Ae(z) is the perturbed dielectric 
profile. All unperturbed states have the same imaginary 
part of their wave vectors (see the inset in Fig. 2) and 
thus their fields have all the same envelope, exponen- 
tially growing outside the slab, with the higher-n states 
oscillating more rapidly, see Fig. [5] (a). In the perturbed 
system, the envelopes are different due the varying Im x v . 
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FIG. 7: Root mean square deviation of the RS electric field 
E„ from its exact value £^ exact ) versus the state number v 
(a) and normalized state number f/n ma x (b), calculated for 
different basis sizes N. The straight magenta line in (b) is a 
cubic law fit. 



Also, as can be seen in Fig.[6f b) , the frequency of the os- 
cillations increases in the perturbed (denser) layer, and 
their amplitudes change at the same time. 

The perturbation theory fully reproduces the EF of the 
RS, both inside and outside the slab. Inside the slab, the 
EF is given by the expansion in Eq. © with the coeffi- 
cients c nu diagonalizing the matrix in Eq. (|7|) . Outside 
the slab, the fields are given by Eq. (0} in which the per- 
turbed wave vectors x v assign the proper exponential 
growth and oscillations of the EF in vacuum, while the 
amplitudes A^j are found by comparing Eqs. (@| and © 
and using the continuity of the EF through the bound- 
aries. To quantify how well the perturbation theory re- 
produces the EF of a RS, we calculate its root mean 
square (RMS) deviation within the system defined by 



A„ = 



C a \£l N \z)-£i cxact) (z)\ 2 dz 



\ f_ a \£i c ™ ct \z)\ 2 dz 



(25) 



The results are shown in Fig. [3 where we have matched 
exact and perturbed RS using the MA and plotted A„ 
for different basis sizes TV. We see that the trend in accu- 
racy with state number and the basis size is the same as 
in Fig.|3](e),(g), and the RMS deviation versus the nor- 
malized state number also shows a universal dependence 
similar to those in Fig.[3](f),(h). However, the EF is in 
general less well reproduced than the wave numbers and 
the power law A v cx (is/N) 3 is observed only in the in- 
terval of 0.05 < A„ < 0.2. 



C. Green's function and transmission 

The Green's function (GF) is an important quantity 
which fully characterizes the response of an optical sys- 
tem, determining its scattering and transmission. For the 
slab with a wide perturbed layer given by Eqs. (|10[) and 
(f22j). the GF G(z, z 1 ; k) which satisfies the equation 

\d 2 z + [e(z) + Ae(z)]k 2 ^G(z : z';k) = 5(z - z') (26) 

and outgoing boundary conditions can be calculated an- 
alytically. Note that when calculating observables, k is 
real as it is given by the vacuum wave number of an exter- 
nal driving field. The GF is calculated using its spectral 
representation 



1.12.13 



£„(z)£„(z') 
2k(k — x v ) 



(27) 



in which the EF £ v {z) and the RS wave numbers yc v 
are calculated numerically via the RSE. For the wave 
numbers x v , we use the extrapolated values Eq. (|20[) . 

In light of the importance of the GF and its further 
usage for calculation of observables, we compare 
the GF calculated by RSE with basis size N and Eq. (|27| . 
to its exact analytic form G'( cxact ) ; again using the RMS 
deviation as given by 



A GF = 



HaS-a G (A °(^')-G( exact 


){z,z>) 


2 dzdz 1 


pa pa 
J —a J —a 


G*( cxact )(z, z') 2 dzdz' 



(28) 

Such a comparison is shown in Fig. [5] for different basis 
sizes N. Increasing the basis size has two effects on the 
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FIG. 8: The root mean square deviation in the GF A GF as a 
function of the wave number of the driving field k, calculated 
via the RSE for different basis sizes TV. 

GF: (i) it improves the GF error at a given k and (ii) 
widens the fc-range of the GF with small error. The latter 
is due to a larger wave-number range of poles in the GF, 
Eq. Ip7)l. being reproduced for large N. 

Both expansions Eqs. ([6]) and (f27|) . for the EF and for 
the GF, are valid only inside the slab or on its borders and 
are not suitable for the vacuum area where the EF of the 
RS grow exponentially. The GF itself is, however, regular 
on the real fc-axis. Moreover, in vacuum, it always has a 
simple analytic form of a plane wave with the amplitude 
that can be deduced from values inside the slab, Eq. (|27p . 
using the continuity of the GF when passing through the 
interfaces. In this way, the GF can be calculated at any 
point of the (z, z') space, inside or outside the slab. 

The delta-function in Eq. ([26]) plays the role of a source 
of plane waves generated at the point z' and propagating 
in both directions, away from the source. The GF then 
has the meaning of the system's response on such a plane- 
wave excitation. This can be used to derive a formula for 
the transmission in terms of the GF. To do this, we place 
the source of strength 2ik just outside the slab at z' = 
—a, in order to produces two plane wave of amplitude 1. 
One of these waves is transmitted trough the slab, and 
just after the slab at point z = a the intensity of the 
EF (which does not change with further increase of z) is 
given by 

T(k) = \2kG(a,-a;k)\ 2 (29) 

and is called transmission. 

We calculate the transmission using Eqs. (|27p and (|2"9")l 
for the GF taken to be either numerical or analyti- 
cal G( cxact ) . This allows us to calculate the absolute error 
in the transmission, |TW — T^ exaoi ^\, which is shown in 
Fig. [51(b). The transmission itself is shown in Fig.^a) 
and has a profile which is fully determined by the pole 
structure of the GF. The RS which contribute in this 
frequency range can be seen in the inset to Fig. 2. We 
see that the error of the transmission has a similar mag- 
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FIG. 9: (a) Light transmission through the slab with a wide- 
layer perturbation Eq. (|22[) . (b) Absolute error in the trans- 
mission calculated using the analytic form of T(k) and nu- 
merical values from the RSE for two different simulations. 

nitudc and scaling with N as the GF itself, as can be 
expected from Eq. (|29|) . 

D. 5-perturbation 

We now move from a wide perturbation to a very nar- 
row and strong one, like a thin metal film on a dielectric. 
Such a perturbation is described by 

Ae(z) = we d S{z - a/2) (30) 

with the delta-scatter strength wed = —0.1a. Physically, 
this perturbation corresponds to a thin layer of the dielec- 
tric constant changed by ed, which is placed at z = a/2 
and has a width w much narrower than the shortest wave- 
length of the resonant modes used in the basis. The di- 
electric profile for the system with the 5-perturbation is 
shown in Fig. [TO] 

As in the case of a wide-layer perturbation considered 
in Section TlV Al we plot and compare in Figs. 11-14 the 
RS wave numbers, calculated exactly and via the RSE 
with and without extrapolation, as well as the parame- 
ters of the power law fit and relative and absolute errors 
which we also need for the quality check of our simula- 
tion and extrapolation. The analytic solutions for the 
(5-perturbation and its matrix elements are given in the 
Appendix [X] 

We see in Fig. [TT] that the extrapolation reduces the 
relative error by 1-2 orders of magnitude. The integral 
strength of the perturbation is much (almost two orders 
of magnitude) weaker than in the case of the wide layer 
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FIG. 10: Dielectric constants of the unperturbed slab e(z) 
and a slab with a 5-perturbation e(z) + Ae(z). The distance 
z is in units of the half width a of the slab. 



considered in Section UV Al However, the convergence is 
much slower in the case of the 5-perturbation. We see 
in Fig.Q2£c),(d) that for large N the power law expo- 
nent is close to a v — —1. This is to be expected as the 
(5-perturbation does not have a finite width. The ma- 
trix elements V nm , though oscillating, have no decrease 
with increasing wave number (or index n) which leads to 
a much stronger mixing of states compared to the wide 
layer perturbation. Indeed, in the wide layer case, states 
with higher indices are less important due to the rapid 
oscillation of their wave functions, so that the matrix el- 
ements scale as V nm cx 1/n (for n ^> m). Using the 
second-order Rayleigh-Schrodinger perturbation theory 
and the explicit form Eqs. (| A9|) and (IA12[) of the matrix 
elements V nm , we can show that the wave number cor- 
rections scale as 1/N and 1/N 3 for the 6- and wide-layer 
perturbations, respectively, in accordance with Figs. 3(d) 
and 12(d). 
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In the case of the 5-perturbation, the absolute errors 
shown in Fig.[T2Tf.h) as functions of the normalized state 
number do not display any universal curves, still for small 
u/N approaching asymptotically a cubic law in the state 
number v (magenta lines). Thus we conclude that in 
this case M v cx v 3 /N [compare with Eq. (|23p]. At larger 
values of u/N this dependence transforms into a linear 
one, M v cx v/N (blue lines). Because of the slow (1/N) 
convergence, the extrapolation gives a huge improvement 
as is clear from Fig.[T3] and demonstrates its necessity in 
the particular case of the 5-perturbation. 

At the same time, the relative extrapolation error is 
predicted within an order of magnitude, as can bee seen 
in Fig.EEl For the majority of RS, i^ cxact) < F v , the 
exact values j7 , j oxact ) being significantly overestimated. 
However, for a large class of solutions it turns out to 
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FIG. 11: As Fig. 2, but for the 5-perturbation shown in 

Fig.rro] 



FIG. 12: As Fig. 3 but for the 5-perturbation shown in Fig. QUI 
The horizontal magenta lines in panels (c) and (d) are a = — 1 
lines. 
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FIG. 13: As Fig. 4, but for the 5-perturbation shown in Fig. 10 
and state number v = 28. 



be highly underestimated. The systematic deviation 
seen Fig. [T?] in estimating the relative extrapolation error 
though F v may be a result of the systematic variation in 
the power law exponent a„ well seen in Fig. ll2( c).(d). 
Hence it is generally advisable when studying conver- 
gence with our method to run simulations with a variety 
of Ni parameters in order to establish over what range 
of Ni the power law is applicable for the given strength 
of perturbation. 

We were also able to simulate a (^-perturbation out- 
side the perturbed slab by taking the unperturbed slab 
to include the position of the delta scatterer and thus 
the perturbation consisting of a superposition of a 5- 
perturbation and a wide layer compensating the differ- 
ence in the dielectric constants between the vacuum and 
the unperturbed slab. In this case we did obtain conver- 
gence of the perturbed wave numbers to the exact solu- 
tion. However, for a c5-perturbation outside of the unper- 
turbed slab or exactly on the border, the simulation does 
not converge to the correct solution. This is to be ex- 
pected since in this case the perturbed RS contain waves 
reflected from the external perturbation, which are waves 
propagating towards the slab. Such incoming waves are 
not part of the basis of unperturbed RS, and thus cannot 
be reproduced by an expansion in this basis. 



E. Microcavity 

To evaluate the RSE in presence of sharp resonances, 
we use a Bragg-mirror microcavity (MC), which con- 
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FIG. 14: As Fig. 5, but for the 5-perturbation shown in 
Fig. 10. 
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FIG. 15: Dielectric profiles of a planar microcavity having 
P — 3 pairs of Bragg mirrors on each side (blue line) and an 
unperturbed dielectric slab (orange line). 



sists of a Fabry-Perot cavity of thickness Lq and refrac- 
tive index jiq surrounded by distributed Bragg reflectors 
(DBRs). The DBRs consist of P pairs of dielectric layers 
with alternating high (tih — 3.0) and low (til = 1.5) re- 
fractive index. In order to a have sharp cavity mode at a 
given wavelength Ac, these alternating layers have to be 
of quarter wavelength optical thickness, and the optical 
thickness of the cavity has to be a multiple of half the 
wavelength. We take Lc — Ac/2. An example of the 
dielectric profile of such a system with P = 3 is shown in 
Fig.[T5l 

The RS of a MC are calculated using the RSE. The 
RS wave vectors and the transmission through the MC 
are shown in Fig. ll6f a).fb). For reference, the unper- 
turbed eigenvalues are also included in Fig.lToTa). The 
unperturbed system taken for the RSE is again a dielec- 
tric slab which dielectric constant e{z) can be seen in 
Fig. [HI Throughout this section the outer boundaries 
of the MC and the unperturbed slab coincide, and we 
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FIG. 16: (a) Wave vectors >c v of the resonant states of a mi- 
crocavity with P = 3 pairs of Bragg mirrors on each side cal- 
culated via RSE with N = 801. (b) Microcavity transmission 
as a function of the normalized wave vector of the incoming 
light; Lc and nc are the cavity thickness and refractive index, 
(c) The difference in the transmission calculated via RSE and 
using the scattering matrix method.— 



choose e s — 5.5 which is between n\ and n 2 H , providing 
good convergence. e s could be further optimized for best 
convergence of the RSE. In order to verify the transmis- 
sion calculated by the RSE, we use the scattering matrix 
method 21 which is a straightforward and precise way of 
calculating the optical properties of a planar system. Fig- 
ure PTBT c) demonstrates a good agreement between the 
two calculations. 

Clearly, there is a one-to-one correspondence between 
the RS wave vectors in Fig. llfif a) and the MC transmis- 
sion in Fig. 1 16( b). Namely, the real part of the wave 
vectors corresponds to the positions of the peaks in the 
transmission while the imaginary part gives their line 
widths. This is well understood in view of the spectral 
representation of the Green's function Eq. (f2~T)) used for 
the calculation of the transmission via Eq. (|29|) . 

One of the modes shown in Fig.fTBTa) is rather isolated 
and has imaginary part much smaller than the others. 
This mode, %c, satisfies the Fabry-Perot resonance con- 
dition Re Kc — n/(L c nc) and is called the cavity mode. 
For the wave vector k of incoming light close to this res- 
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FIG. 17: The FWHM (a) and the position of the cavity mode 
(b) calculated analytically and via the RSE for different num- 
ber of pairs P of Bragg mirrors on each side of the microcavity. 
N is the basis size used in the RSE. Where possible, extrap- 
olated wave numbers have been used. Crossed rectangles for 
TV = 51 indicate states which are rejected by the SC. 



onance condition, k « ir/(Lcnc), the Greens' function 
Eq. (|27|) is dominated by a single term corresponding to 
this narrow mode. As a consequence, there is a sharp 
peak in the center of a wide stop-band seen in the trans- 
mission in Fig. irTJT b'). For sufficiently large P an ana- 
lytic approximation for its full width at half maximum 
(FWHM) is known,^ 2 - 



Afc 



4n c 



2P 



1 



L c + 



A, 



n L n H 



(31) 



2 nc{n H -n L ) 



which we use to compare with the RSE calculation. With 
the refractive index of the external material n cx t = 1 
and using Ac = 2Lc and nc = tin, Eq. (|3ip reduces 
to Ak = 4(n# — riL)(nL /nn) 2P / {Lcn 3 H ) . Comparison 
of the above formula with the RSE result for the cavity 
mode is given in Fig. 1171 for different number of Bragg- 
mirror pairs P and for different basis size N in the RSE. 
Figure [T3 demonstrates that RSE is capable of giving 
both the correct width and location of sharp resonances 
in the transmission profile, if a large enough basis is used, 
in spite of there being no sharp resonances in the basis. 
As the basis size is enlarged, the width and the peak lo- 
cation of the cavity mode converge to the analytic values. 
The fact that for a fixed N the cavity mode position and 
the width are predicted worse for larger P is explained by 
our choice of the unperturbed slab which always has ex- 
actly the same thickness as the Bragg-mirror MC. With 
the number of Bragg-mirrors increasing, the field inside 
the MC oscillates more rapidly (also shifting the cavity 
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mode towards higher frequencies) that requires a larger 
number of RS to be taken into account in order to pro- 
duce results on the same level of accuracy. We have ver- 
ified (not shown) that the errors become independent of 
P, if one and the same constant width of the unperturbed 
slab is used for different values of P. 



V. SUMMARY 

The resonant state expansion has been implemented 
and validated in planar open optical systems reducible 
to effective one-dimensional systems. A reliable method 
of calculation of resonant states, and in particular their 
wave numbers, electric fields, as well as the Green's func- 
tion and the transmission of such systems, has been de- 
veloped and demonstrated^ It includes estimation of the 
accuracy and convergency of calculations and in particu- 
lar extrapolation of the eigen-wavevectors towards their 
exact values which are generally not available. Particu- 
lar examples which illustrate the general method and the 
developed algorithm include a dielectric slab with wide- 
layer and ^-perturbations as well as an optical microcav- 
ity having different number of Bragg mirrors. In these ex- 
amples, a comparison with exact solutions has been made 
in order to verify the approach. In all three systems the 
resonant states and the transmission are reproducible to 
any required accuracy by the resonant state expansion. 
The wave vectors of resonant states are the most essential 
part of the calculation as they most strongly affect the 
optical properties of the system through the poles of the 
Green's function. The extrapolation of the wave vectors 
using the power law in the basis size, which has been de- 
veloped and demonstrated, significantly improve the ac- 
curacy of calculations, by one or two orders of magnitude. 
Application of the method to two and three-dimensional 
systems will be reported in future works. 
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all having the same imaginary part. 



(A4) 



2. Resonant states of a slab perturbed by a wide 
dielectric layer 

The exact solutions of the wave equation Eq. ([3]) for 
the system with the perturbation given by Eq. (|2"2"j) and 
outgoing boundary conditions has the form 

A ve -ixuz 5 z < -a, 

.(exact) M= J ^e'^ 2 +ae-'^ z , -a<z<b, 
» Dve*^*"* + E^e-V^* , b < z < a , 
H„e l>c " z , z > a , 

(A5) 

where e p = e s + Ae, and b — a/2. The coefficients in 
Eq. (IA5j) are found from the continuity of the electric 
field and its derivative and the normalization condition 
Eq. (|24p . The complex- valued RS wave numbers x v are 
found by solving a secular equation following from the 
boundary conditions: 



hf{k)g(k) - 1 = 



/3-7 
07-1 



09(h)- if (k) 



(A6) 



where 



Appendix A: Analytic calculation of resonant states 
and perturbation matrices 



(A7) 



1. Resonant states of the unperturbed slab 

Solving the wave equation Eq. (|3]) with Ae(z) = 
and the profile of the dielectric constant e(z) given by 
Eq. (fTUj) . the electric field of RS n, normalized according 
to Eq. ([5]), takes the form 



and the functions f(k) and g(k) arc defined as 

f(k) = e -2iv^*<«+i0 , g{k) = e -HsfcKa.-b) ( A g) 



We solve Eq. (|A6I) using the Newton-Raphson method to 
findfc = ^ cxact) . 
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3. Matrix elements of the wide-layer perturbation takes the form 



Using Eq. (JSJ and basis functions Eq. (|A1[) we calculate 
V nm for the wide-layer perturbation Eq. (j2"2")l to be 



1 



e s Ma^/T s 



+ (-z) n - m r 1 (k n -k m ,z) 
+ {-i)- n+m r 1 {-k n + k m ,z) 

+(-z)- n - m V (-k n ~k m ,z) 



(A9) 



for n m and 




■1) 



n [77(2fc n ,z) + 77(-2/c„,z)]° 



4iau/e 



for n = to, where ^(fc, 2;) = e l ^^ kz /k. 



(A10) 



4. Resonant states of a slab perturbed by a delta 
scatterer 

In the case of a (^-perturbation Ae(z) = wed5{z — b) 
with \b\ < a, the secular equation for the RS wave vectors 



It is also solved numerically with the help of the Newton- 
Raphson method to find k = x^ act \ 



5. Matrix elements of the 5-perturbation 



Using Eq. (JSJ and basis functions Eq. (IA1|) we calculate 
V nm for the 5-perturbation to be 



V nm = we d E n (a/2)E m (a/2) 



(A12) 
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